module sphere_geometry_mod

  use flogger
  use const_mod
  use math_mod

  implicit none

  private

  public euler_formula
  public lonlat2xyz
  public xyz2lonlat
  public rotate
  public rotate_back
  public great_circle
  public spherical_area

  interface lonlat2xyz
    module procedure lonlat2xyz_r8
    module procedure lonlat2xyz_r16
  end interface lonlat2xyz

  interface xyz2lonlat
    module procedure xyz2lonlat_r8
    module procedure xyz2lonlat_r16
  end interface xyz2lonlat

  interface rotate
    module procedure rotate_r8
    module procedure rotate_r16
  end interface rotate

  interface rotate_back
    module procedure rotate_back_r8
    module procedure rotate_back_r16
  end interface rotate_back

  interface great_circle
    module procedure great_circle_r8
    module procedure great_circle_r16
  end interface

  interface spherical_area
    module procedure spherical_area_3_r8
    module procedure spherical_area_3_r16
    module procedure spherical_area_4_r8
    module procedure spherical_area_4_r16
  end interface spherical_area

contains

  integer function euler_formula(nc, nv, ne) result(res)

    integer, intent(in), optional :: nc
    integer, intent(in), optional :: nv
    integer, intent(in), optional :: ne

    if (present(nc) .and. present(nv)) then
      res = nc + nv - 2
    else if (present(nc) .and. present(ne)) then
      res = ne - nc + 2
    else if (present(nv) .and. present(ne)) then
      res = ne - nv + 2
    end if

  end function euler_formula

  subroutine lonlat2xyz_r8(R, lon, lat, x, y, z)

    real(8), intent(in ) :: R
    real(8), intent(in ) :: lon, lat
    real(8), intent(out) :: x, y, z

    real(8) cos_lat

    cos_lat = cos(lat)
    x = R * cos_lat * cos(lon)
    y = R * cos_lat * sin(lon)
    z = R * sin(lat)

  end subroutine lonlat2xyz_r8

  subroutine lonlat2xyz_r16(R, lon, lat, x, y, z)

    real(8 ), intent(in ) :: R
    real(16), intent(in ) :: lon, lat
    real(16), intent(out) :: x, y, z

    real(16) cos_lat

    cos_lat = cos(lat)
    x = R * cos_lat * cos(lon)
    y = R * cos_lat * sin(lon)
    z = R * sin(lat)

  end subroutine lonlat2xyz_r16

  subroutine xyz2lonlat_r8(R, x, y, z, lon, lat)

    real(8), intent(in ) :: R
    real(8), intent(in ) :: x, y, z
    real(8), intent(out) :: lon, lat

    lon = atan2(y, x)
    lat = asin(z / R)

    if (lon < 0.0d0) lon = lon + pi2

  end subroutine xyz2lonlat_r8

  subroutine xyz2lonlat_r16(R, lon, lat, x, y, z)

    real(8 ), intent(in ) :: R
    real(16), intent(out) :: lon, lat
    real(16), intent(in)  :: x, y, z

    lon = atan2(y, x)
    lat = asin(z / R)

    if (lon < 0.0d0) lon = lon + pi2

  end subroutine xyz2lonlat_r16

  ! ************************************************************************** !
  ! Rotation transform                                                         !
  ! Purpose:                                                                   !
  !   Calculate the rotating transformation and its inverse of the original    !
  !   coordinate system (lon_o,lat_o) to the rotated one (lon_r, lat_r) with   !
  !   the north pole (lon_p,lat_p) defined at the original coordinate system.  !
  ! ************************************************************************** !

  subroutine rotate_r8(lon_p, lat_p, lon_o, lat_o, lon_r, lat_r)

    real(8), intent(in) :: lon_p, lat_p ! Rotated pole coordinate
    real(8), intent(in) :: lon_o, lat_o ! Original coordinate
    real(8), intent(out), optional :: lon_r, lat_r ! Rotated coordinate

    real(8) tmp1, tmp2, tmp3, dlon

    dlon = lon_o - lon_p
    if (present(lon_r)) then
        tmp1 = cos(lat_o) * sin(dlon)
        tmp2 = cos(lat_o) * sin(lat_p) * cos(dlon) - cos(lat_p) * sin(lat_o)
        lon_r = atan2(tmp1, tmp2)
        if (lon_r < 0.0d0) lon_r = pi2 + lon_r
    end if
    if (present(lat_r)) then
        tmp1 = sin(lat_o) * sin(lat_p)
        tmp2 = cos(lat_o) * cos(lat_p) * cos(dlon)
        tmp3 = tmp1 + tmp2
        tmp3 = min(max(tmp3, -1.0d0), 1.0d0)
        lat_r = asin(tmp3)
    end if

  end subroutine rotate_r8

  subroutine rotate_r16(lon_p, lat_p, lon_o, lat_o, lon_r, lat_r)

    real(16), intent(in) :: lon_p, lat_p ! Rotated pole coordinate
    real(16), intent(in) :: lon_o, lat_o ! Original coordinate
    real(16), intent(out), optional :: lon_r, lat_r ! Rotated coordinate

    real(16) tmp1, tmp2, tmp3, dlon

    dlon = lon_o - lon_p
    if (present(lon_r)) then
        tmp1 = cos(lat_o) * sin(dlon)
        tmp2 = cos(lat_o) * sin(lat_p) * cos(dlon) - cos(lat_p) * sin(lat_o)
        lon_r = atan2(tmp1, tmp2)
        if (lon_r < 0.0d0) lon_r = pi2 + lon_r
    end if
    if (present(lat_r)) then
        tmp1 = sin(lat_o) * sin(lat_p)
        tmp2 = cos(lat_o) * cos(lat_p) * cos(dlon)
        tmp3 = tmp1 + tmp2
        tmp3 = min(max(tmp3, -1.0d0), 1.0d0)
        lat_r = asin(tmp3)
    end if

  end subroutine rotate_r16

  subroutine rotate_back_r8(lon_p, lat_p, lon_o, lat_o, lon_r, lat_r)

    real(8), intent(in)  :: lon_p, lat_p ! Rotated pole coordinate
    real(8), intent(out) :: lon_o, lat_o ! Original coordinate
    real(8), intent(in)  :: lon_r, lat_r ! Rotated coordinate

    real(8) sin_lon_r, cos_lon_r, sin_lat_r, cos_lat_r, sin_lat_p, cos_lat_p
    real(8) tmp1, tmp2, tmp3

    sin_lon_r = sin(lon_r)
    cos_lon_r = cos(lon_r)
    sin_lat_r = sin(lat_r)
    cos_lat_r = cos(lat_r)
    sin_lat_p = sin(lat_p)
    cos_lat_p = cos(lat_p)

    tmp1 = cos_lat_r * sin_lon_r
    tmp2 = sin_lat_r * cos_lat_p + cos_lat_r * cos_lon_r * sin_lat_p
    ! This trick is due to the inaccuracy of trigonometry calculation.
    if (abs(tmp2) < eps) tmp2 = 0.0d0
    lon_o = atan2(tmp1, tmp2)
    lon_o = lon_p + lon_o
    if (lon_o > pi2) lon_o = lon_o - pi2
    tmp1 = sin_lat_r * sin_lat_p
    tmp2 = cos_lat_r * cos_lat_p * cos_lon_r
    tmp3 = tmp1 - tmp2
    tmp3 = min(max(tmp3, -1.0d0), 1.0d0)
    lat_o = asin(tmp3)

end subroutine rotate_back_r8

  subroutine rotate_back_r16(lon_p, lat_p, lon_o, lat_o, lon_r, lat_r)

      real(16), intent(in)  :: lon_p, lat_p ! Rotated pole coordinate
      real(16), intent(out) :: lon_o, lat_o ! Original coordinate
      real(16), intent(in)  :: lon_r, lat_r ! Rotated coordinate

      real(16) sin_lon_r, cos_lon_r, sin_lat_r, cos_lat_r, sin_lat_p, cos_lat_p
      real(16) tmp1, tmp2, tmp3

      sin_lon_r = sin(lon_r)
      cos_lon_r = cos(lon_r)
      sin_lat_r = sin(lat_r)
      cos_lat_r = cos(lat_r)
      sin_lat_p = sin(lat_p)
      cos_lat_p = cos(lat_p)

      tmp1 = cos_lat_r * sin_lon_r
      tmp2 = sin_lat_r * cos_lat_p + cos_lat_r * cos_lon_r * sin_lat_p
      ! This trick is due to the inaccuracy of trigonometry calculation.
      if (abs(tmp2) < eps) tmp2 = 0.0d0
      lon_o = atan2(tmp1, tmp2)
      lon_o = lon_p + lon_o
      if (lon_o > pi2) lon_o = lon_o - pi2
      tmp1 = sin_lat_r * sin_lat_p
      tmp2 = cos_lat_r * cos_lat_p * cos_lon_r
      tmp3 = tmp1 - tmp2
      tmp3 = min(max(tmp3, -1.0d0), 1.0d0)
      lat_o = asin(tmp3)

  end subroutine rotate_back_r16

  pure real(8) function great_circle_r8(R, lon1, lat1, lon2, lat2) result(res)

    real(8), intent(in) :: R
    real(8), intent(in) :: lon1
    real(8), intent(in) :: lat1
    real(8), intent(in) :: lon2
    real(8), intent(in) :: lat2

    res = R * acos(min(1.0d0, max(-1.0d0, sin(lat1) * sin(lat2) + cos(lat1) * cos(lat2) * cos(lon1 - lon2))))

  end function great_circle_r8

  pure real(16) function great_circle_r16(R, lon1, lat1, lon2, lat2) result(res)

    real(16), intent(in) :: R
    real(16), intent(in) :: lon1
    real(16), intent(in) :: lat1
    real(16), intent(in) :: lon2
    real(16), intent(in) :: lat2

    res = R * acos(min(1.0d0, max(-1.0d0, sin(lat1) * sin(lat2) + cos(lat1) * cos(lat2) * cos(lon1 - lon2))))

  end function great_circle_r16

  real(8) function spherical_area_3_r8(R, p1, p2, p3) result(res)

    real(8), intent(in) :: R
    real(8), intent(in) :: p1(:)
    real(8), intent(in) :: p2(:)
    real(8), intent(in) :: p3(:)

    real(8), dimension(3) :: p31, p12, p23
    real(8) a312, a123, a231

    p31 = cross_product(p3, p1); p31 = p31 / norm2(p31)
    p12 = cross_product(p1, p2); p12 = p12 / norm2(p12)
    p23 = cross_product(p2, p3); p23 = p23 / norm2(p23)
    a312 = acos(-dot_product(p31, p12))
    a123 = acos(-dot_product(p12, p23))
    a231 = acos(-dot_product(p23, p31))
    res = R**2 * (a312 + a123 + a231 - pi)

  end function spherical_area_3_r8

  real(8) function spherical_area_4_r8(R, p1, p2, p3, p4) result(res)

    real(8), intent(in) :: R
    real(8), intent(in) :: p1(:)
    real(8), intent(in) :: p2(:)
    real(8), intent(in) :: p3(:)
    real(8), intent(in) :: p4(:)

    real(8), dimension(3) :: p41, p12, p23, p34
    real(8) a412, a123, a234, a341

    p41 = cross_product(p4, p1); p41 = p41 / norm2(p41)
    p12 = cross_product(p1, p2); p12 = p12 / norm2(p12)
    p23 = cross_product(p2, p3); p23 = p23 / norm2(p23)
    p34 = cross_product(p3, p4); p34 = p34 / norm2(p34)
    a412 = acos(-dot_product(p41, p12))
    a123 = acos(-dot_product(p12, p23))
    a234 = acos(-dot_product(p23, p34))
    a341 = acos(-dot_product(p34, p41))
    res = R**2 * (a412 + a123 + a234 + a341 - pi2)

  end function spherical_area_4_r8

  real(16) function spherical_area_3_r16(R, p1, p2, p3) result(res)

    real(16), intent(in) :: R
    real(16), intent(in) :: p1(:)
    real(16), intent(in) :: p2(:)
    real(16), intent(in) :: p3(:)

    real(16), dimension(3) :: p31, p12, p23
    real(16) a312, a123, a231

    p31 = cross_product(p3, p1); p31 = p31 / norm2(p31)
    p12 = cross_product(p1, p2); p12 = p12 / norm2(p12)
    p23 = cross_product(p2, p3); p23 = p23 / norm2(p23)
    a312 = acos(-dot_product(p31, p12))
    a123 = acos(-dot_product(p12, p23))
    a231 = acos(-dot_product(p23, p31))
    res = R**2 * (a312 + a123 + a231 - pi)

  end function spherical_area_3_r16

  real(16) function spherical_area_4_r16(R, p1, p2, p3, p4) result(res)

    real(16), intent(in) :: R
    real(16), intent(in) :: p1(:)
    real(16), intent(in) :: p2(:)
    real(16), intent(in) :: p3(:)
    real(16), intent(in) :: p4(:)

    real(16), dimension(3) :: p41, p12, p23, p34
    real(16) a412, a123, a234, a341

    p41 = cross_product(p4, p1); p41 = p41 / norm2(p41)
    p12 = cross_product(p1, p2); p12 = p12 / norm2(p12)
    p23 = cross_product(p2, p3); p23 = p23 / norm2(p23)
    p34 = cross_product(p3, p4); p34 = p34 / norm2(p34)
    a412 = acos(-dot_product(p41, p12))
    a123 = acos(-dot_product(p12, p23))
    a234 = acos(-dot_product(p23, p34))
    a341 = acos(-dot_product(p34, p41))
    res = R**2 * (a412 + a123 + a234 + a341 - pi2)

  end function spherical_area_4_r16

end module sphere_geometry_mod
